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SnCTIO.N 1 


ItvTRODUCTION 


The concept of using a solar sail as a means of transportation 
vfithin the solar system was first discussed seriously in the United 
States in the 1950's.^ Recently, the poss'bility of using a solar 
sail has been given serious consideration. The earliest contender for 
an actual mission v.'ould be a rendezvous with Halley's comet in 1986, 
with a launch in 1982. After that a solar sail could be used for a 
variety of solar system missions including, for example, a Mars sample 
return. Hov;ever, recent NASA decisions indicate that the solar sail 
v/ill not be developed for near term missions. Some possible missions 
would include a planetary escape or capture segment of the trajectory 
in addition to a heliocentric trajectory segment. The thrust acceler- 
ation level for sail designs in the near term are low, on the order 
of 1 mm/s , and so escape or capture requires a spiral phase from a 
low orbit to escape or from initial capture to a lov^er orbit. Typi- 
cally there would be several dozen or even hundreds of orbits. This 
system is .similar to other low thrust systems such as electric pro- 
pulsion. The principal subjects of this report are the spiral phase 
of a planetary escape trajectory and orbit to orbit transfer. A 
computer program has been developed to calculate solar sail planetary 
trajectories to a near escape point. It can also be used to calculate 
spirals to a lower orbit and for trajectories about Mercury, Venus, 
and Mars as well as Earth. 

a 

A solar sail consists of a large area of reflective material. 

Typically, a very thin plast .e is coated with a reflective metallic 

surface, Sunlight reflecting off the surface imparts a force which 

is , approximately normal to the surface. The pressure exerted on a 

2 

perfectly reflecting sail normal to the sun at 1 A.U. is 0.9 N/m . 

By varying the angle of the surface v;ith the photon direction, the 
force direction can be changed, although there will never He a 
component of force toward tJie sun. If the sail has a large enough 
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aroa, and the total vehicle is light enough/ then reasonable acceler- 
ations will result. In heliocentric space it is possible to thrust 
BO that the sail jnoves away from the sun by having a eomponent of the 
reflective force in the direction of motion. If the force has a 
component in the direction opposite to the motion, the total angular 
momentum will be decreased and the sail vehicle will fall in toward 
the Dun, When the sail is in orbit about a planet, the sail direction 
can be varied so that the thrust almost always has a component in the 
direction of motion and so will increase the semimajor axis or energy 
even when the vehicle is moving toward the sun in its plunetocentrie 
orbit, 'Typically, planetary orbits have periods of only a few hours 
or days and, therefore, the necessary sail angle v;ith respect to the 
sun direction to obtain desired force directions may be changing 
rapidly. In practice there are constraints on how fast these sail 
angles can be changed. These constraints are (jonsidered only in a 
limited fashion in this report. 

The Jet Propulsion Laboratory has considered two sail models as 
possible interplanetary propulsion systems: a square sail and the 

heliogyro. More recently a decision was made to concentrate study 
on the heiiogyro model. This report stresses the square sail model 
which wa.<f implemented in the computer code. A heliogyu model is 
discussed which is similar to the square sail model, but which was 
not implemented. 

The square sail is a large (about 850m X 850m) , nearly flat 
sheet with attendant structure, one side reflective and one side dark. 
The heliogyro consists of a central hub with 12 blades, 6 each in 
two parallel planes, each blade about 6 km long and several meters 
wide. The heliogyro would spin with a period of about three minutes. 
Rigidity is provided by the centrifugal force. Variations in thrust 
direction as Well as torques for attitude changes are accomplished by 
varying the blade pitch, that is, the angle about the longitudinal 
axis of the blade. 

The purpose of this study was to develop a ptJgram to produce 
optimal planeooentric solar sail trajectories and to apply that ” 
program to a performance analysis. Earlier papers ahve considered 
nonoptimal planetary trajectories. Sands^ considered two-dimension-^ 
al escape maneuver for a trajectory plane containing the sun-planet 
line for the idealized flat sail, perfectly reflecting specularly on 
both sides. The sail is assumed to rotate about its axis at half the 
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rate of revolution about the planot. Pimple*' conuider'id the throe-' 

dimensional naae of a sail vehicle in orbit about a planet with the 

orbital plane normal to the aun-planot line, also Cor tlio idealised 

Bail. The nail angle is chosen so that the component of thrust in 

the direction of motion iri maximized. Quasi-circular orbits result 

until near escape . Neither Sands nor Pimple included solar motion 

Time optimal helicentric trajectories have been calculated for the 

idealized sail, Zhukov and Lobedev** considered planar traj^-Jtories , 

Sauer^ generalized this work to three-dimensional trajectories. 
n 

MaclJeiL conceived of ti^o heliogyro design and with others produced 
additional studies . ® ^ 

The study reported here built upon a computer program (SECKSPOT) 

which had been developed to produce electric propulsion geocentric 
n 1 !? 

orbit transfers. ' The program uses; (1) the method of averag- 
ing in order to reduce the amount of computer time heeded to 
calculate many trajectories, each of which includes many orbits about 
a planet, and (2) equinoctial orbital dlements which are nonsingular 
for zero eccentricity and inclination. Since averaging is valid 
only if the thrust to local weight ratio is small, the scheme^ cannot 
by itself yield trajectories to escape energy (infinite semimajor 
axis) . The scheme can be used to calculato trajectories to a large 
semimajor axis (small magnitude but negative energy) . The effect of 
shadowing by the planot may be included as v/ell as oblateness for 
Earth trajectories. The solar sail program (called SUNSPOT 'or SUN- 
Sail Program for Optimal Trajectories) does not include planetshine 
or drag effects which may be important at low altitudes. Although 
attitude constraints may be important for some trajectories, they 
were not included explicitly in the optimization procedure. 

Kryloff-Bogoliubof f averaging of both the state and costate is 
used. The averaged rates of change of the mean values of the state 
and costate are found by numerical quadrature. The differential 
equations for the mean state and costate may then be integrated in 
large time steps (typically days). The method of averaging has 
been used extensively in recent years, Edelbaum'*' ' “ has used 
averaging to calculate analytic solutions for special cases of optimal 
low thrust trajectories, and others have used averaging when con- 
sidering effects such as oblateness, third body perturbations and 
non-optiinal thrusting. ' Jasper'^'^ utilizes equinoctial orbital 
elements and averaging in lov; thrust optimization work. The effect 
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of oblateneas is included by analytically adding its associated rate 
of change of the moan state and costate to that duo to thrust. The 
effects of shadowing are calculated by assuming that there iis no solar 
radiation pressure when the sail is in shadow. 

The overall trajectory is optimized by a shooting method. 

Initial values of the unspecified states and costaUes are chosen at the 
initial time. An optimum low thrust trajectory is then generated by 
integrating the state and costato until the final time. This 
will generate an optimal trajectory to the wrong terminal state. A 
sensitivity matrix is then generated by varying the initial conditions 
and running a set of neighboring trajectories. A Newton iteration on 
the initial conditions is then used to drive the terminal errors to 
within specified bounds. The final converged trajectory is a minimum 
time trajectdry ( except when a penalty function is added to the cost 
to prevent subterran' trajectories, in which case nearly minimum 
time trajectories result) . 

A computer program has been developed to calculate planetocentric 
solar sail trajectories. The analysis and code can be used in a per- 
formance analysis. A limited set of cases are discussed in the result 
section of this report. A paper based on the material in this report 
has been presented."^ Related work for trajectories beginning at 
large semimaj or axis and continuing to escape energy performed by Green 
is reported in Ref. 22. The results of the work reported here and of 
the work by Green include the first production of optimal solar sail 
planetocentric trajectories. 

In the next section some of the general techniques used are dis- 
cussed. The following section contains the main analytical contribu- j 
tion of this study. A number of trajectories are given in the result 
section, constituting a partial performance analysis. A preliminary 
heliogyro model is discussed in Appendix A. Results from previous 
efforts which are needed for the sail program are given in Appendices 
B - D. As part of the solar sail effort, preliminary results yielded 
non-optimal n«Jar escape trajectories using a maximum energy strategy. 
These results were reported to JPL during the course of the work, but 
are summarized in Appendix B. 


SECTION 2 


GENERAL TECHNIQUE 


2.1 Introduction 

Three areas of general technique are discussed in this section. 

One is the method oi’ averaginr*, essential to the running of a program 
which generates many trajectories in a reasonable amount of computer 
time. Second is discussed the method of generation of the final 
trajectory, h time optimal trajectory is desired, A state and costato 
formulation is used which results in a two-point boundary-value problem 
which can be solved by a Newton iteration procedure. Finally some 
comments on numeri.cal techniques are made, 

2.2 Averaging 


A great savings in computer time can be effected by considering 
a first approximation to the state and costate. Short period varia- 
tions in the state and costato are eliminated by the averaging tech- 
nique. When low thrust propulsion is utilized and the other 
perturbations to the inverse square motion are small and when the state 
includes the five slowly varying orbital elements which indicate 
the size, shape and orientation of an orbit and possibly other slowly 
varying quantities, then averaging may be used. The orbital element 
indicating the position of the spacecraft in the orbit is eliminated 
by the averaging process . - 

The averaged Hamiltonian can be defined as 


H = i 2 H at 

^ ffl 


{ 2 . 1 ) 


where H is the unaveraged Hamiltonian and T is the orbital period. 
When calculating this integral the state and costate are held fixed. 
The motion of the spacecraft is assumed to vary in a manner described 
by Kepler's equation over the averaging integration. The approximate 
state and costate satisfy the canonical equations, 
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( 2 . 2 ) 






where the overbar in.licafciis the approximate quantities. 


(2.3) 


In what follows the averaging integral for oblateness (J^) is 
solved analytically; otherwise a numerical quadra t\sre formula is used. 
The differential equations can then be solved numerically using a 
time step which is much larger than but unrelated to the number of 
orbital revolutions. 


2,3 Solution of Boundary Value Problem 

The method used to generate a low thrust trajectory is to develop 
the Hamiltonian, to calculate a control (the thrust direction) and to 
write the canonical equations for the state and costa te. The initial 
state is specified. Applicat.’on of transversality conditions for the 
time optimal problem yield aduitional specification b on the state and 
costate. Thus a two-point boundary-value problem results which must 
be solved to obtain the requisite trajectory. When these equations 
are solved, an extremal trajectory will result which is usually locally 
optimal. No attempt is made to investigate generalized Jacobi"type 
conditions to establish local sufficiency. Also, in common with other 
nonlinear problems, there may be more than one extremal meeting the 
same end conditions. The very difficult question of global optimality 
is not considered. 

The single trajectory generation portion of the code is coupled 
with a Newton iterator to solve the two-point boundary-value problem. 
The unknown initial conditions and value of the final time are iterated 
on in order to meet the final conditions which are functions of the 
final state and costate. The partial derivative matrix of final condi- 
tions with respect to the initial costate-is obtained numericully 
by calculating neighboring trajectories to a nominal. 

The Newton method works by first guessing values for the iteration 
parameters, call them x and t^^, and then running a nominal trajeobory 
which will yield final conditions which in general are not equal to 
the desired final conditions, y^. Rev1.L.,-d values- for X, t^ may then 
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be obtainud by caXculafcing a scnoitivifcy matrix or partial derivative 
matrix, A, which is generated by varying slightly, one at a time, 
each of the iteration parameters, x, and running a new, neighboring, 
trajectory. Differencing the resulting values of the final conditions 
With the nominal values yields a &y for each in addition can 

bo oulculutod analytically except pQr!^aps for ^ which can bo ^ 
approximated numerically by varying tg slightly and evaluating 
the corresponding H, differencing this with the nominal H and dividing by 
At^, Then A is an approximation for the partial of y, with respect" 
to X, tfl. 


( 2.0 


A revised estimate of the iteration parameters can then be obtained by 
the formula 


NEW 


— ■ - 
X 

•=£ 


- '‘"bz-ia’ 


(2.5i 


A new nominal trajectory can then be generated and the procedure 
continued until the final conditions are met to within some tolerance. 
In the event that the nev; x, t^ do not yield a reduction in the norm 
of the final condition errors, the change in (x, t^) is reduced in 
magnitude by factors of 2, Also there is an option of using a modi- 
fied Newton-Raphson procedure wherein the A matrix is not always 
recalculated at each iteration by running neighboring tvajectories , 
but instead a new A may be approKimated using the old A and the 
values of the changes in x, t^. 

2 . 4 Numerlral Methods 

A Newton-Raphson iterator is used which calculates the sensitivity 
matrix by running neighboring trajectories by changing slightly the 
initial valvies of the iteration parameters, one at a time. The size 
of the change in the iteration variables is chosen by the user and 
can affect the accuracy of the matrix, A Modified Nev'ton-Raphson 
iterator used basically the same technique but many of the iterations 
make use of a modified sensitivity matrix rather than calculating a 
new one by running neighboring trajectories at each itv'sration. 
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The low thifUBt differuntial equations are integrated using a 
fourth order Runge-Kutta method. The time stop is selected by the user. 
Cutting the size of the time stop can increase the accuracy of the 
trajectory but rapidly increase run time. 

Numerical averaging utilizes a Gaussian quadrature. The number 
of points sampled on an orbit can largely be determined by the user. 
Again, more points increase accuracy at the expense of run time. 


SECTION 3 ^ 

ANALYSIS V 


3 . 1 Sail Force Model 

A square sail force model was implemented in the computer code and 
this section stresses that model. Limited heliogyro data was supplied 
by dPL and preliminary modeling effort is discussed in A peridix A, This 
effort was not included in the computer code. A simple .model of the 
heliogyro would be similar to that given in this section. j;' 

The acceleration exerted by the photon pressure of sunlight on a 
flat, perfectly x'cflectlng surface, whose normal has an angle, a, with 
the sun direction, has a direction along the normal to the baek^of the 
sail and a magnitude given by 


“c 2 " ■ 

a„ = — g- cos a (3.1) 

R 


where a_, the characteristic acceleration, is the acceleration at 1 A.U. 

Q 

caused by the photons reflecting off a sail whose surface is normal to 
the sun-lino; R is the distance from the sun in AiU.'s. If R is the unit 
vector from sun to space vehxcle and n the unxt normal to the sail back> 
the acceleration can be written 




“c 

- (n-^R), n 


(3.2) 


since cos a = n R. This expression represents the model that has been 
used in much of the previously reported work {e.g., Ref. 3, 4, 5, 6). 

The square sail model which was under consideration was not 
' perfectly flat, nor was the reflection purely specular. The framework 
produced a kite-like structure. The sail was bowed somewhat with an 


apt'X tho I'cntcr. Tttcic wax nom** abnurption <ind somr vmiasion from both 
front and back. To moto iccuratoly moil«l thusp u( foots, ai. uxpa xiun .n 
cox t*, whore 6 is the .in<|le k>etwipti the force vector and tlie B'in>llne, 
was fit to more precise data for the ijeometry and ri fli*ctlve charactcr- 

1 3 

istics r.f the sail model by thi‘ Jet Propulsion Laboratory.* The follow- 
ln<j relationship was produced] 


** A 

rtp • (Cj ^ 0*2 cos 20 ♦ Cj cos 46)u (3.3) 

K 

A 

wh#*re u is unit vector in the force direction, not necuBuurily p'.allel 

A 

with the nutmal to nail nurfAce. The thrust dirertion u could >e written 
as a function of two angles, the cone angle, 0, wher«* cos 0 • H***u and the 
cloc); angle, which in the angle between projections of a roforenco direc- 
tion and the thrust direction onto a plane normal to tlie sun-lino (see 
Figure 3-1). The expression in Fq. (3.3) could also bu expressed as a 



Figure 3-1. Sail geometry (in the plane of the sun-line 
and force vector) . 
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2 4 

function of cok 0 and con fj , The parancturs c^,C 2 «c^ are constants which 

depuna on the shape and reflective qualities of the sail. Note that if 

2 

c^ C 2 ■ 0.5 end Cj ■ 0, then the exfression reduces to cos i , which 
correspond:, to the idealized square sail. The relationship between a 
l>ody>fixcd cooidinate frame and 5 may be obtained but is not necessary 
for our purposes. 

One model produced by JPL had ccefficlents given by ■ 0.367, 

Cj ■ 0.643, and c^ ■ -0.010. The locus of the tip of th^ normalized 
force vector for these coefficients and also for the idealized sail 
coefficients is plotted in Figure 3-2. A line drawn from the origin at 
an angle 0 with the horizontal axis and terminating at the curve will 



Figure 3-2. Locus cf tip of force vector as cone angle varies. 
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have a lenyth repteaenUing the force magnitude. The maximum force magni- 
tude is one, which occurs when the force vector is parallel to the sun- 
line, At 0 e 62.585® the force magnitude is zero for the nonidealize^ 
sail. This la the cutoff angle. Por 6 greater than this angle, the 
equation dictates a negative force magnitude, j 

The force model of Eq, (3.3) is only an approximation and its be- 
havior for 6 approaching 180® is not realistic. If c^ ^ 0, then for some 
0 < 90® I b(0) “ 0 and for larger 0, b is negative implying a force toward 
the sun, which is impossible. The 0 at which b(D)= 0 was considered a 
cutoff point and any 0 which resulted from the mathematical optimization 
greater than this cutoff was reset to the cutoff value. The cutoff value 
may be obtained by setting b(0) = o and solving for & analytically, as a 
function of Cj^, c^, and c^. 

In order to calculate the cutoff angle ii. is convenient to trans- 
form Eq. (3,3) (o an expression in cos 6 and cos^O, i.e. 

ii " ' ■ 

|i b(0) “ + )C2 cos^O + cos^O (3.4) 

Prom, trigonometric identities: 

ki = Cl - C2 + C3 

k2 “ 2C2 - 8C3 (3.5) 

^3 = 


For b(0) - 0 


- k- 




cos 0 = 


4kik3 


2k. 


(3.6) 


The coefficient kj was always positive and much larger than k^ end k^ 
so that the positive square root is the correct root. Thus, the 
cutoff value of 0 is given by 



cos 



A; - <kj^>=3 

2 k, ’ 


(3.7) 
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whoro Oj, ia between 0“ and 90® (in the computer code 0 and, therefore, 

0_ waa assumed to lie betv/een 90® and 180®; thus, the siqn on the outer 
c 

square root was negative) . If « 0, the cut-off angle is simply 
obtained, 

^ ^ The Equations of Motion in Equinoctial Orbital Elements 

) I ri 

A variation of parameter formulation using equinoctial orbital 
elements was used for the equations of motion. By using equinoctial 
orbital elements the singularities that occur for zero eccentricity or 
inclinations of zero or ninety degrees when using classical orbital 
elements are avoided. (For inclinations near 180®, retrograde equinoctial 
orbital elements can be used, although we will not consider that case 
in this report.) 

The direct equinoctial orbital elements are defined in terms of 
the classical orbital elements by the formulas 

a “ a 

i'l 

h = e sin (o» + fi) 

■ : " -= Q cos (m + 9.) — - (3.8) 

p « tan(^) sin 9 

q “ tan (^-) cosn 

where a is the semimajor axis, e is the eccentricity, i is the inclin- 
ation, SI is the longitude of the ascending node, m is the argument of 
pericenter. A sixth parameter, F, called the eccentric longitude, indi- 
cates ,the position in an orbit. It is given in terms of the classical 
variables by 


P = E + n + m (3.9) 

Where E is the eccentric anomaly. The variable F will be eliminated 
by the^averaging process. Further details about equinoctial elements 
are given in Refs. 12, 17 and 24, 

The inverse relationships are defined by 

= a 

= (h^ + 
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e 


x.vsr i 


A 




(3.10) 


n ■ tan (p/q) 

b> - tan”^(h/)t) - tAn"^(p/q) 


T)ie equinoctial coordinate frame is defined by t)>e basis vectors 
f, w, which are given below with respect to an inertial coordinate 
frame. 

fl - 


1 ♦ p‘ + q 


1 + p‘ + q 


2 2 
1 ♦ p - q 


(3.11) 


w - 2 J - 2 q 

1 T p + q 2 2 

|l - P - q 


This coordinate frame is illustrated in Figure 3-3 where w is normal to 
the orbital plane. 



ORBITAL 

PLAr.E 

i .\ / 




UNIT 

SPHERE 


Figure 3-3. The equinoctial coordinate frame. 
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The mean longitude is defined by 

A - ;t ♦ w •» n 


(3.12) 


The eccentric longitude, T, v.as defined in Zq. (3.9). Kepler's equation 
in terms of X and r is then given by 




X 

- F - 

k sin F ‘f h cos F 

(3.13) 

and 

velocity 

are given by 





r ■ 

Xii ♦ Vji 

(3.14) 




• 

r “ 

^li + ^i2 

(3.15) 


m 

a( (1 

- h^g) 

cos F + hkB sin F - kl 

(3.16) 


m 

aid 

- 

sin F ♦ hkB cos F - hj 

(3.17) 

^1 

- 

na! 

r 

[hkS cos 

F - (1 - h^B) sin F) 

(3.18) 


- 

nai 

r 

1(1 - 

?) cos F - hkB sin FJ 

(3.19) 




n 

• 

(3.20) 





a 




r 

A 

- 1 - 

k cos F - h sin F 

(3.21) 


M is the gravitational const/.nt and 


S - 

1 + ,4 - h^ - )c^ 


(3.22) 


omciNAL MCI? to 

p(xm ql'auty 
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The unaveruged variation of parasnefcors formula for an object in 
an inverse square, gravitational field perturbed by a force u is 

,i S B M(2,F) u (3.23) 

where Z represents the five equinoctial orbital elements (a, h, kf p, q) . 
We define the 5x3 matrix 

32 , 

M(Z,F) B -7 '• (3.24) 

3r- 

The elements of this matrix are given in Table B~1 of Appendix B, 

3« 3 The Optimization Problem 

It is desired to calculate minimum time trajectories for escape, 
capture, and orbit to orbit transfer. The orbital element and averaging 
techniques did not permit trajectories to escape energy when eccentri- 
city goes to one and the thrust to weight ratio becomes too large. For 
planetary escape the final state condition is that energy be zero. This 
is equivalent to infinite semimajor .axis, Because of the limitations of 
the technique, trajectories to a subescape condition are considered. 

This is defined as a large, but finite semimajor axis, or equivalently 
as a small magnitude, negative energy. The initial orbit is assxiroed to 
be given. For orbit transfer the initial orbit is given, and either 
the final orbit is completely specified or else three orbital elements 
(a, e, i) are specified. The capture problem is considered a special 
case of the orbit transfer problem, sincd a zero energy initial orbit 
cannot be assumed. Thus a quasi-capture trajectory problem assumes that 
the initial orbit is given and that the final orbit is at some lower 
energy requiring a spiral trajectory downward. 

A calculus of variations or maximum principal approach is used. 

The initial time and state are given; some subset of the final state is 
given. The state equations are given by Eq. (3.23) with the force given 
by Bq. (3.3), In summary 

a 

Z = -§• b(e) M(Z,F) u (3.25) 
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a » 


whoro 


- fi 0' 


b((t) ej + Cg C0B2Q + Cg coh40 


(3.26) 


GosO « R u 


(3.27) 


anti Gj^, 02/ arc copstantB. The distance to the sun is a function of 

time and the direction R depends on the state when using an equinoctial 

coordinate frame. The cliaractoristic acceleration, a„, is a given 

/> ^ 

constant. The control is the force direction u. 


The unaveraged Hamiltonian is given by 


fl», m A 

1 Z ” b ( 0 ) X Mu 
- - 


(3,28) 


Applying tho method of averaging, define the averaged Hamiltonian as 


r ^ 

— I H(f, X, F(t), t, u) dt 

o' ■ ~ 


(3.29) 






where T^ is the orbital period. It is convenient to perform the 
integfation with respect to tho eccentric longitude, F, Then 


ir 

H = ™ y* II (Z, X, F, t, u) (ll) & 


(3.30) 


where 


5F “ 2^" (1 - X cosp, - FT sinF) 


(3.31) 


For" convenience define 


= ^ 3T 


(3.32) 
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The overbuy .indicates that the quantities are held constant during the 
averaging interval. The t dependence is time dependence which does not 
depend on the spacecraft location such as the planet - sun position. 

Necessary conditions ^for an optimal trajectory are that the 
Eulor-Lagrange aquations be satisfied; 


where 







(3.33) 


(3.34) 



Z (Z, F, G, t) 


(3.35) 


3H 

3Z, 


<r A ~2T 

“ 3Z, 


M + b^ 

R /sz^ DZj^ 


u s 


(3.36) 


3£ 

3Z 


1 

-tv 


0 

~ sinP 
- cosF 
0 
0 


(3.37) 


3M * 

- — is given in Appendix B. 
3Z 


-0 
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I'tom ligo. (3.26) amt ^3, 27) 


9b „ 9b 90 


az 


n 


where 


9 b 


■g^ - 2 g 2 sin 2S - Ac^ sin ;-40 


and 


H 

93 


1 G'^ 


3R 


sin6 


az 


BuU since sin29 “ 2sin0 cosG and sin46 « 4sin0 cos0 cos20 then 


3 

~ “ (;4c, COS0 + 16c, cos6 cos20) G"^ — 

az \ 2 ^ “ az 


V 


If R is given in equinoctial coordinates where R = (X , Y , Z_) 

h’ 5 S S " 


the nonzero partials are given by 


ax 


5 

ap 


iq'i ^ Z ) 

1 + ® ^ 


ax. 


2pY 


8q 


5 

2 5" 


1 + p + q 


ax. 


2qX^ 


3p 


1 + p + q 


Ik 


9Y, 


oq.. 


1 + p + q 


r— r? + ”s> 


6mG<NAi<. to 

( OB fom quaijt^ 


ill! 

sl'l 

I 
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I 


(3.38) 


(3.39) 


(3.40) 


(3.41) 

then 


(3.42) 


(3.43) 


(3.44) 


(3.45) 






1 + p'' + 


O.'l?) 


The control must maximize the namilhonian. That portion of tho 
Unaveraged Hamiltonian containing the control is 


H' « b(G) u 


(3.40) 


Mow b is considered a function of u rather than 0. It is convenient to 
identify the primer vector «ib 


iv " - 


(3.49) 


and define the pri.'tjer vector cone angle $ where 


A 

cosS « R Xy 


(3.50) 


The Hamiltonian is then 


H' « b(u) Xy u 


(3.51) 


The control u can be written as a function of two angles, the cone 
angle 0 and a cloclc angle, tj;. Than 


G “ (cost^ ej^ + sin'J» Gg) sinO + cos8 e^ 


(3.52) 


v/here the unit vectors are given by 


B.3 ^ S 


(3.53) 


—2 “ — 3 ^ y~ 


(3.54) 


El “ £2 ^ ^3 


(3.55) 


f 


It 


where V ia a roCorenco vector, in t)iis coordinate frame the primer vector 
can bo given in terma of B and a clock angle Then 

« (coflijijj + ainijjp Sg) ain& h- cooB o^ (3.56) 

The relevant portion of tlie Hamiltonian ia then 

U' “ b( 0 ) (coBi|) cosijjp H- oihij' sint)»p) ainO ainp + cosO coaB (3.S7) 

Maximizing H' with reapect to 6 and ^ is equivalent to maximizing H' 
with respect to u. The clock angle for u must be equal to the clock 
angle of the primer vector; i|< » »}<p. Thus R, and u are in the 
same plane and from geometry; 



sin(B - e) 
'iln'g 


R + 


ainO ? 

SInF iv 


(3.58) 


(as in Ref. 6 ). 


The II' reduces to 

H' » (Cj^ + C 2 COB 20 + C 2 COB 46 ) (cosO cosB sinO sinB) (3.59) 


Thi‘< can be maximized with respect to 0 by setting its derivative with 
respect to 0 to zero. 

a (- 2 G 2 sin26 - 40^ sin 46) (cos6 cos3 + sin9 sing) 


+ (c^ + C2 cob2D + c., cos 40) (- sinO aiii3 + oosO sin5) 

= 0 (3.60) 

Now for thei idealized flat square sail = Cg 0.5 and “ 0 so that 

b(6) ■= cos^O. Then w = 0 can bo reduced to a quadratic in tanO whose 

d 0 

solution is 


tanO 




3cos3 + /8sxn 3 4 
4sin3 


9cos^B 


(3.61) 


original pagh; i& 

OF POOR QUALITY; 


21 


whuru note that P can vary between 0 and lUO®, but 6 varies between 0 
and 90“ . This is similar to nof. G. in the general case, tliu oppression 
could be reduced to a polynomial, but it is easier to solvo the eguatioa 
numerically using a Newton method. If 

II « b(0) d(0) 

(3.62) 

where 


d(0) cosO cosg + slnO sin? 

(3.63) 

lot 


f(0) u » b’(0) d(0) + b(0) d'(0) 

(3.64) 

The prime indicates the derivative with respect to 0. An 
for 0 is given by the solution for the idealized sail, Eg. 
iterate 

initial guess 
(3.61). Then 

'’n+1 “ 

(3.65) 

Until' “ Ojj < c, some small positive number. Also needed is 

fi ~ >= 2b'd' + b" d + b‘d" 

UD 

(3.66) 

where 


b" =s - 4 c2 cos 20 - cos40 

(3.67) 

and 

• 


d" = -d 

(3.68) 

A 

Thus the control u is determined. 


Figure 3-4 shows the variation of 0 with S for the idealized sail 
and for one set of coefficients; c^^ = 0.367, “ 0.643, c^ “ -0.010 

which corresponds to one, JPE supplied, heliogyro model. Note the 
cut off angle of 0 » 62,6“ which occurs in a region where the model is 
not realistic. In the computer program, whenever the optimization 





PRIMER VECTOR CONE ANGLE (d»fl> 

F'igura 3-4. Thrust vector cone angle versus 
primer vector cone angle. 


called for a 9 greater than this cutoff value. 6 was set to the cut- 
off value as Indicated in the figure. The cutoff angle was given by 
Eq. (3.7). 

The state and costate lifferential equations have been given and 
the optimal control derived. It re.mains to tpecify the transversality 
conditions. For transfer to a subescajre point, the final energy, or 
equivalently, the final seminajor axir. is specified. 


E 


m 

I 


n u 

2 a 


(3.69) 


This is equivalent to specifying the final semimajor axis since the 
mass, m, is constant. Transversality conditions then require that 
the adjoints to the remaining ocate variables be zero and that the 
Hamiltonian bo equal to one at the final time. 
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0 


0 

0 




- 1 


(3.70) 


The adjoint* to the classical orbital eler.ertts, e, i, fi, and u should 
also be zero at the final time. 

In the case of orbit transfer, if the urbit is completely speci- 
fied (a, e, i, n, u or equivalently a, h, k, p, q) at the final time, 
then their adjoints are free and the only transversal.ity condition is 
that the final Hamiltonian be one. If three orbital elements (the 
semimajor axis, eccentricity, and inclination) are specified at the 
final time, then the Hamiltonian must be one and and must be zero. 
In terns of equinoctial elements and their adjoints the following final 
state conditions are specifiedt 


a “ 

/pW' 

The transversality conditions aret 


a (t j) 

’ e(tj) 

- tan(i(t^)/2) 


hX^ - kXj, 


- 0 


P\j - ^^p “ “ 0 


(3.71) 


(3.72) 


The state and costate differential equations and the initial and 
final conditions yield a two-point boundary-value problem. This problem 
can be solved with en iterative method. A Newton method was used as 
in Ref. 12. The differential equations were integrated numerically with 
a Runge-Kutta method. At each function evaluation for the Runge-Ku'eta 
method a numerical quadrature was used to average the equations (3.33) 
and (3.34). The control was found using Eq. (3.58) and by a Newton search 
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for the cono iimjlo. At3difcional complexity ia added by including a pen-^ 
ality function i tshadowing and oblatenenB. 

3. 4 Perrcentor Penalty Function u " , - 

Planetoceiitric solar sail,,-trajectoi ;ies typically build up large 
eccentricity rapidly. This can' be particularly undesirnblo if the ini- 
tial altitude is low. Optimization may actually yield trajectories that 
intersect the planet's surface. Also, large eccentricity orbits often 
require rapid changes in thrust direction which, in practice, may be 

!l 

impossible to implement. One way of reducing this problem is to append 
a penalty function to the cost (which has been flight time), Ohe 
possible penalty is the integral over the flight time of the inverse 
pericenter squared. Then the cost is given by 


n ? i 
2 J 2 


(3,73) 


This penalizes low altitudes in general and especially at small semi- 
major axis. To some extent it also penalizes larger eccentricities. 

The constanj:, n, must be picked to obtain a de^.irable weighting between 
the penalty and the flight time. , Experience with a few example cases 
indicated .that this method could be used tc prevent altitudes which 
were too sltiall with a very small increase in flight time (less than 5%) . 
Other penalty functions might be used, for example, one which is much 
stvOeper near some specific altitude, thus acting as an inequality 
cjjnstraint (see Ref. 25) . 

Jj 'The Hamiltonian and the, costate equations must be modified. The 
pericenter is given in equinoctial elements by 


a(l - 


+ k^) 


(3.74) 


2 . , 2 , 


a"(l - /h" + k") 


(3.75) 
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tlien the Hamiltonian is given by 


H « -=L 


(i;7G) 


The averaged Hamiltonian is 


H B 




L(Z) dt 


(3.77) 


^ind the oostate equation is 


I 


(3.78) 


The first integral , -in Eq.iij(3.77) is the same as if there were no penalty 
function. The second intii|gral does not include a time dependence in ''7 
the integrand and, so it is' just equal to L(£). Thus, the only thing 
different in Eg. (3.78) is'^^the addition of ^ , In particular '■,1,, 


'V- 8L 

\ ,1 3 n 


(3.79) 


3L' 

IK " 


1 - /h^ + 


(3.80) 


1 - Ai.^ + 


(3,81) 


^ = ill 
Dp 3q 


(3.82) 


Transversalit;y conditions still require that the final Hamiltonian be 
equal to one, but in this case at t.; 


3,5 OblattiiieBs 


In tho previous sections we have considered only perturbations to 
tile inverse square motion caused by thrusting, in this section, the 
effect of oblatoness (J 2 ) is considered, Oblateness equations are in- 
eluded here for convenience from Ref. 12. Oblateness is included only 
for earth-centered trajectories. 


Tho single-averaged perturbing potential due to has been cal- 
culated in terms of equinoctial coordinates in Ref. 14 and is repeated 
in Appendix C. R^ is the equatorial radius of the Earth and J 2 is set 
to .001827. These formulas enter the averaged Hamiltonian as coeffi- 
cients of the costate (outside the integral since the averaging effect 
has already been accounted for) . 

If z indicates the perturbation due to thrust as given in Eq. (3.33) 

>> 

then the Hamiltonian is given by 


fpJU mj. 

- X^Z_ + X'^Z, 
— — j 2 — ~a 


(3„84) 


The state equation is 


Z = Z_ + z^ 


(3,85) 


The costate equation is 


= 


3Z^ 

3H _ _ tT ~^2 


3Z 


3Z 


a 

-I 

-7T 


X- s + X Z ^ 
3Z 3Z 


dF (3,86) 


The partials indicated by 3Z_ /3Z in the above expression are given in 

— j 2 — 

Appendix C. 


3.6 The Shadow Effect 

For solar sail missions, there is no thrust while the spacecraft 
is in a planet's shadow. The entry and exit angles are needed in order 
to perform the averaging integral. in calculating these angles the 
following assumptions are made. The shadow is cylindrical; the planet 
revolves around the sun in an elliptical orbit; and over one spacecraft 
revolution, the sun's direction is fixed. Pertinent equations for the 
calculation of the entry and exit eccentric longitudes are summarized 
in Appendix D taken from Ref. 12. 
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Let F, refer to the eccentric longitude at entry to the shadow 
and at exit. That part of the Hamiltonian proportional to thrust 
is then 



(3.87) 


It 


(3.88) 


(3.89) 


The calculation of — is discussed i,n Appendix D. 

dZ " 

3.7 Planetary Data and Coordinate Frames 

Earlier versions of the code for electric propulsion trajectories 
assumed Earth orbital missions with spacecraft orbital elements refer- 
enced to an equatorial coordinate frarne,._ The srilar sail program has 
been generalized to include trajectorita 'about the four inner solar 
system planets: Mercury, Venus, Earth, and Mars. The gravitational 

constant for each planet is given in Table 3-1, Also gwen is the 
planet’s radius. This is used in shadow calculations and also the 
"internal units" of the code use planetarv radii. Oblateness, is 

assumed zero for all planets except Earth. The characteristic accelera- 
tion of a solar sail is defined as the maximum acceleration at 1 A.U. 
Thus, the planet’s distance from the sun is needed as well as the sun's 
direction. An equatorial coordinate frame is used only for Earth tra- 
jectories, The obliquity angle, the angle between the Bartli's equator 
and the .ecliptic, is thus heeded. For all the planets an ecliptic frame 
may, be used. In this frame the X- and Y-a>:es are in the ecliptic with 


0 




ORIGINAL PAGE 
OF POOR quality 


. li 

-■ ■ (l 

the X-aNifj towuctl UIig veirnal equinox ditecfclon and the X-^axis is normal 
to the ocliistic in n “nox'therly" direction. Anotlier coordinate frame 
option ia the "planethrv" frame which is referenced to the planet's 
orbital plane, the K-axis normal to the orbital plane, the X-axis pointed 
toward perihelion. 


Table 3-1. Planetary data. 



Mercury 

Venus 

Earth 

Mars 

Semimajor axis, a (A.U.) 

0.3Q7099 

0.723322 

1.0 

1.523691 

Eccentricity, e 

0.205627 

0.006793 

0.016726 

0.093368“ 

Inclination, i (“) 

7.00399 

3.39423 

0, 

1.84991 

Longitude of the ascending 
node, n(“) 

47.05714 

76.31972 

0. 

49.24903 

Mean longitude of perihelion 
u {") 

76.83309 

131.00331 ■ 

102.25253 

335.32269 

Mean orbital motion, 
n ("/day) 

4.092339 

1.602131 

0.905609 

0.524033 

Mean longitude at epoch, e 

222.62165 

174.29431 

100.15815 

258.76729 

Obliquity angle (") 

— 

-- — 

23.45 


•^2 

— 

— 

0.0010827 


Radius (km) 

2435.0 

6052.0 

6378.16 

3393.4 

Gravitational constant, 
U (km3s“2) 

22181,6 

i24860,l 

// 

398601.2 

42828.4 


I’i 


Data for Epoi^h 1960 Jan 1.5 E.T. (J.D. 2436935.0) from Ref. 26 except J 2 , 
radius and ji from Ref. 27, 


In the ecliptic coordinate frame tlie unit vector pointing fx'om the 
planet toward tile sun is 


“S 



cosH 

-cosi sinf) 

sini sinn 


-cos (v+uj) 


sinn 

cosi cosH 

-sini cosH 


-sin (v+m) 


0 

sini 

ccsi 


0 






.. .. 


(3.90) 
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Note thut for Barth, 1 and n are zero so^ tliat the matrix is the identity 
matrix. The argument of perihelion u is obtained from Tablo 3-1 since 

uJ « w - n T. (3.|l) 

0 

The angle v is the true anomaly. It can bo approximated by (Ref. 2G) 


V a M + (2e -^) sinM + ^ sin2M ^ e^ sin3M 


{3v-9.2) 

■A 


where M is the mean anomaly and 


M e nt + 


(3.93) 


The time^ T is measured from the epoch and is the mean anomaly at 
epoch; 


Mq e - 01 


(3.94) 


The mean orbii)al motion, n, is given in the table with other needed 
constants. ~ " 

The distance to the sun in A.U. 's is 


l^sl 


a(l - e Q 
1“ + e cosv 


(3.95) 


The difference between the sun-spacecraft distance and the sun-planet 
distance is assumed negligible. ■ 

- -" For an Earth-equatorial frame 


-cos (tu+v) 
-cosO sin(w+v) 
“SinO cos (w+v) 


(3.96) 


where o is the obliquity of the ecliptic-. 

For a "planetary” frame i, 0, and w in Eq. (3.90) are set to zero, 





In tliQ uquinqetial coojfdinate framq 


( 3 * 1 ?) 

where the equinoctial basis vectors were given in Eq. (3,11). Note 
that the R vector used in earlier equations in this section is the 
vector from the sun to the planet rather than from the planet to the 
sun . ' 


^ ^ ^ A |« A 

Bg “ il< H) £g 
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SKCTTON 4 


NUMKHICAL R^ISULTS 


A number of example runs were performed. A complete perfoimanco 
analynis was not done because of the contlnuln<j charyu in status of the 
Holar anil mission planning. First ttie square sail was emphasized. 

Then th9 hcliogyro was picked <is the principal solar sail design. Fi- 
nally, the solar sail lost in competition with ion drive as a contender 
for near future low thrusl miss.'.ons. A complete heliogyro model was not 
implemented. Kesulv.s in this section are for the square sail model 
discussed in Section 3. A limited number of Car'.h-centered trajectories 
will be discussed hero. Because of decreased interest, further runs 
were not (>erformcd. Ko trajectories about the other inner planets arc 
included since the only cases that were run w<;re short test cases to 
verifv the coding. A test case also verified the use of the program 
to calculate trajertozies which spiral down from a high near-capture 
orbit to a lower orbit. This is a special case of orbit transfer. 

A nu..iber of cases uning the square sail apprex i ir.etion with helio- 

gyro coefficients ■ 0.367; C 2 ■ 0.643; Cj ■ -0.010) for transfer 

to a Bubescape point will be discussed. These runs used a 10 day time 

step and four 4-point gaussian quadratures for the averaging integral. 

A subcscape point with a semimajer axis of 100,000 or 200,000 km was 

2 2 

used (this is equivalent to a Cj of -3.99 or -1.99 km /s , respectively). 

A number of factors determine flight time. The initial orbit is, of 

course, of critical importance. An initial orbit with a semimajor 

axis of 21378 km, eccentricity of 0.655 was suggested by JPL personnel. 

This initial orbit would be representative of a tug-launched sail 
2 2 

(Cj ■ -18.6 km /s ). The longitude of the ascendzng node (0) and the 
argument of perigee (w) can be varied. The inel i-'-^tion (i) can be 
varied but differing launch energies would be required if t launch 
day is fixed. A numl>er of open loop trajectories w»*re run to get a 
feeling for the i, II, and u that would yield the lower flight times. 

Then tra jf?ctories wore optimized for a few cases of specific i, R, and 

i-teCCOlNG PAa^ Bl^NK NCTT FlL«t0 
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y {for conveHienoQ/ those angloo wero spocified with respect to an 
Gcjliptic frame). A launch of March 21, 1904 was assumed. Although on 
ollipbieal Earth orbit about the sun was modeled, the 1/R effect of 
solar pressure was not included in the exumploa (a constant 1 A.U, dis- 
tance was assumed) . Earlier examples indicated that including oblate- 0 
ness has a very small effect on flight times and therefore this effect 
was not included in the examples shown here. Shadowing was. included 
but thfe Initial orbits were chosen to avoid the shadowed regions. 

The results for seven examples are shown in Table 4-1, cases 1, 

3, and 4 show the effect of varying churacteristlc acceleration from 

0. 6 to 1.0 mm/s . Those cases are also plotted in Fig. 4-1. case 2 
has a final semiraajor axis of «00,000 )?r rather than 100,000 km as in 

'icase 1. An* additional 20 days is needed. At 100,000 km the orbital 

\i, , . ■ 

[period is 3.6 days and at 200,000 km it is 10,3 days, Escape would 
iprobably occur in less than 30 more days, about two more revolutionsV. 
Apooenter for the final' orbit was at 347,000 km where the thrust/weight 

V '*' ' r, -i .0 r. 2 ‘ ■' 

ratio was nearly 0.1 for the 0.6 «un/s acceleration, so that the aver- 
aging technique has limited validity j-there. 

Case 5 is similar to Case 1 except that the initial orbit is cir- 
cular rather than elliptical. The initial C, for both cases is -18.6 
km /s . The transfer time to a of -3.-99-ls about 4 days longer for 

the initially circular orbit case. This behavior is due to the fact 
that the initial orientation of the elliptical orbit was chosen so » 
that the sail was moving away from the sun near apocenter. At this 

r . I ‘i 

point, the sail was moving more slowly and therefore, received more 
energy. Thus, for the first 80 days of Case 1, the eccentricity 
actually decreases to below 0.2. 

Cases 6 and 7 have initial orbits which are normal to the ecliptic 
and chosen so that as the Earth moves around the sun, | the orbits ^will 
come closer to being normal to the sunline. Slightly increased trans- 
fer times resulted compared to the 45® inclination cases. For Case 6, 
the eccentricity became smaller and at the final time was equal to 
.004. The initial costates are shown for all the runs.-;. 

Case 1 will be illustrated in more detail: by a series of figures. 
Figures 4-2 to 4-6 are plots of the classical orbital elements (a, c, 

1, .7, and a) . since the initial time is March 21, the x-axis of the 
coordinate frame is pointing toward the sun at the initial time. 
Eccentricity and inclination decrease for the first 80 days. The line 

- i I ’ .. i\ 
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TABLE 4-1 


TiransfGr to a Subescape Point 


CASE 

1 

2 

- 3 

4 

5 

6 

7 

Initial a (km) 

21378 

21378 

21378 

21378 

21378 

21378 

21378 

e 

0.GS5 

0.655 

0.655 

0.655 

0.0 

0.655 

0.655 

i 

45° 

45° 

0 

in 

45° 

45° 

90° 

90° 

n 

150° 

150° 

150“ 

150° 

150° 


0° 

u 

180° 

180“ 

180° 

180° 

180° 

-90° 

-90° 

Characteristic 
Acceleration (mm/s'^) 

O.fi 

0.6 

0.8 

1.0 

0.6 

0.6 

0.6 

Final a (km) 

100000 

200000 

100000 

100000 

100000 

100000 

200000 

Flight Tima (days) 

lie. 5 

147.4 

88.1 

70.9 

120.3 

122.8 

153. 3 

Costate A ^ 

3389 

3750 

2579 

2097 

3444 

3753 

3653 

Xr 

25S 

84 

-112 

-249 

-801 

12-75 ■ 

1616 


-453 

-464 

-366 

-347 

1365 

66 

66 


-869 

-2614 

-966 

-972 

-603 

763 ' 

804 


1840 

920 

1312 ' 

977 

1059 

-2 

-1 















of nodes changes slightly but the argument of perigee increases from 
-180° to 35°. Figure 4-7 shows a plot of C^. The dotted line at 
“ 0 is the escape energy. Although shadowing was included, the 
trajectory never intersected the shadow. Fairly large changes in sail 
direction were required in a short time for early orbits. For example, 
on the first orbit a change in the .sail direction of over 40° in 20 
minutes was called for. This occurs in the region where the sail is 
near to edge -on with the sunline and the thrust produced is very small. 

In practice a slower non-optiinal change in direction would probably have 
only a small effect on performance. The cone angle for the primor vector 
and the cone and clock angle for the thrust vector are shown in Figs. 

4-8 to 4-10 for 3 orbits including the initial orbit, an orbit at 50 
days, and the final orbit. 

Figures 4-11 to 4-20 show the history of the equinoctial orbital 
elements and their adjoints for Case 1. The initial cosbate is shown 
in Table 4-1. When oblateness was included for the same initial condi- 
tions, the flight time was 116.42 days rather than 116.52 days. The 
corresponding costahe is {3355, 227, -387,-698, 1877), A run with a 
five day time-step was not appreciably different t.han the ten-day case. 
The curves in Pigs. 4-21 to 4-23 have not been smoothed? the ten day 
time step produces corners in the computer pldt. Figures 4-21 to 4-23 
Illustrates the chpnqe- in period", pericenher and anocenter. , ' , 

Convergence characteristics tor the cases of transfer to a sub- 
escape point were good. Typically an initial guess of = 1000 with 
the other adjoints zero resulted in convergence (to within 3 or 4 
significant figures in the desired semimajor axis) in 4 or 5 iterations. 
This requires about 1 minute of CPU time on Draper Lab's Amdahl 470 
(a -fairly fast computer) at a cost of about $20. 

An idealized sail was assumed for an orbit raising case (c^ = C 2 
= 0.5, - 0 ) . A launch time of March 21, 1977 for an initial circu- 

lar orbit with a = 7078 km, i = 28.3°, R = 0°, m s 0 ° in an equatorial 
coordinate frame and with a characteristic acceleration of 1.0 mm/s 
was specified. Shadowing was not included. The desired final orbit 
was equatorial geosynchronous. For this case the minimum transfer 
time was 195.8 days. The eccentricity increased to over 0.83 at 120 
days before decreasing to zero. Unfortunately, the pericenter de- 
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inibially, ryacUing 4060 km afc XIO days, Iji orddi.' bo proven b 
1 / bh?i lowering of purlconter, a penalby funcbiohrwas added bo tlje coab a£i 
disoussed in Secbiqn 3.^ Afber some brial and error adjuabmenb of bhe 
weighbing faobpr bebweon biine and bhe penalby funcbion, a brajucbory 
vms tjenerabed v/ibh a branafer bime of 200.7 daya with a minimum peri~ 
center of 7650 km and a maximum' eccentriei by of 0.72, Thus bhe brans'- 

// _ ■ f. 

far bime penalty waa leas than. ^3? . I'igures 4-?.4. to 4-27 show the aeini- 
major axiOi eccentricity, inciination and pericenbor liishories for- the 
case wibh bhe penalty function plus the paricenter history for the": . 
minimum time case. Although the trajectory was not changed -very much, 
the costate histories were changed considerably. For example, the 
initial costate for the minimum time transfer was (5842, 3070, 3351, I 
-545S, -4979). With the penalby function the initial costate iyas 
(9204, -372, 6533, -7903, ~8536) . It was necessary to perform a series 
of runs with increasing weighting^facUor in order to obtain^i convergence . 
The penalty function might be necessary for low altitude runia with 
shadowing since currently the rhadow code is ndt valid for tr-djectorioa- 
■ which intersect the planet's surface. - ’ 

' ?, These results are meant to illustrate the capabilities of the pro- 
ilgram andl of course, do not" represent" an extensive performance study. " 
i’Typical flight times are indicated, however, and show that the sail 
can reasonably be used for Karbh escape and for orbit transfer. 
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F.igure 4-24 Semimajor axis history 

Orbit transfer with penalty function 


61 












PER I CENTER (KMJ (XID 


fi 


H 


f > 

*trt 


r\i 



•1j.ra 


> 1- 

40. 00 


^ ^ 

60.00 l^l].00 


ll>0.L'C 

TIME 


«) 0.00 
UJRTSJ 




— H 

PUO.CO 
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SECTION* 5 


CONCLUDlN'n REMARKS 


A method for producing optimal solar sail planotocuntric trajec- 
tories has been devised and implemented in a computer program called 
SUNSPOT. A square . ail force model is assumed alt(iough it is similar 
to a simple heliogyro model. Ability to produce orbit to a subescapo 
point and orbit to orbit transfer exists. The method of averaging 
allows rapid computation of trajectories, but since the thrust to 
weight ratio is assumed small, it does not allow trajectories to escape 
energy (infinite semimajor axis) . The subescape point is defined as a 
largo but finite seminajor axis. Typically one or two more revolut*'»ns 
would bo required for escape and a suboptimal scheme can give easy and 
good estimates of the remaining flight time. In order to prevent inter- 
section of a planet's surface, a penalty function may bo included in 
the cost. Trajectories about the four inner planets can bo generated. 

The program is based on an earlier electric propulsion geocentric 
transfer program. As with the earlier worlc, shadowing, oblateness, and 
solar motion may be included. Equinoctial orbit elements are used to 
avoid singularities for zero eccentricity and inclination. 

The optimization requires the solution of a two-point boundary- 
valuu problem. Several trajectories must be calculated to obtain con- 
vergence to the required boundary condi-ions. Each trajectory is cal- 
culated by integrating a first approximation to the state consisting 
of the orbital elements and the costate using a time step of several 
days. Averaging integrals are performed using a gaussian quadrature. 
Typically 16 points are used. In order to obtain convergence in a few 
iterations, a fairly good initial guess is required, especially for 
orbit-to-orbit transfer. Because th< ^olar sail always has a component 
of force away from the sun, eccentricity tends to change quickly. 
Therefore, attainment of a particular final eccentricity can bo diffi- 
cult. Transfer to a subescape point did not present convergence 
dif f icultics . 
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/ 

A tuiiubes of ejiampl0 tr/fij^akorieB v?erG pt'oduoed rtlUhovtyh a thorough 
porfomancG analyBits v/an nijit done because of changittg coindJ fcions in bho 
sail Mission planning, tnibial conditions on tuany of the cases we.ro 
selected to take advantage of geoisietry to avoid the shadow and decrease 
flight time. Examples ate for Earth orbital trajectories. Plight 
times are reasonablo for a low thrust vehicle. 

There are several possible extensioxs of this ’'>»ork. One is the 
inclusion of a more accurate heliogyro sail model. Planetshine effects 
may be important at low altitudes and for Mercury trajectories as indi- 
cated by some preliminary effort for this study. Attitude constraints, 
includinc rate and ra^he change constraints, are important. Such con- 
straints v/ere not included in the optimization although resulting tra- 
jectories can be inspected for violation of the constraints. Typically 
the constraints are /violated during a segment of the orbit v/hen the 
thrust is small combared to the rest of tho orbit and so probably a 
suboptimal pointing vjould have a small effect. There is a tendency 
for eccentricity tjo become large so that intersections of the planet’s 
surface occur. Currently, the shadow computations are not valid if 
this happens, cau/sing ah abort. A penalty function was used to pre- 
vent surface int^irsection for some examples, and there was only a small 
loss of performa/nce Obtaining convergence can be difficult, laspeclally 
for orbit raisi^ig when a small final eccentricity is desired-. For a 
complete escape/ trajectory a precision trajectory integrati-in segment 
could be connected to the initial averaged section with the appropriate 
transition. 
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APPENDIX A 


PPELIMINAHY IIEtilOGYPO FORCE MODELING 


Accurate force modolintj of the heliogyro ia vei'y difficult. The 
force depends on the angular momentum vector, pitch' angle , rotation 
angle, deformation etc. In an effort ho produce a more simple model 
JPL obtained approximate thrust magnitude versus thrust vector cone 
angle values for various sun cone angles (i.e. the angle betv/een the 
angular momentum vector or spin axis and the Bun-Tline) . Oniy data 
for zero collective pitch was obtained. Thus the sun-line, spin axis 
and thrust direction are assumed to lie in a plane. This data is 
reproduced in 'pable A-1 and this case is illustrated in Fig. A-1. 


TABLE A-1 

Values of Normalized Force Versus Force Cone Angle 

r. 

(for two sun cone angles and for 0° collective pitch; 
taken from JPL data) 


Cyellle Pitch 
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r IN SUN SPIN PIANI- 
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SPIN 

AXIS 

Mqure A-l lialioqyro Geometry 

The spin vector direction can chanqe only very slowly. As a 

first cut the dnta for zero sun cone anqle was used to try to obtain 

values for coefficients uslnq the square sail model. Since the data 

2 

diverges considerably from a cos 0 curve, that form did not yield a 
qood fit. A form usinq cos40 anvi cosSa did yield a very qood fit. 

In particular lot 

- 0. 333+0. 77O9COS40-0. 1042 cos80 (A.l) 

.9012 

Then Table A-2 showa a comparison of values from this formula with 
data from Table A-2. Note the cut-off anqle of 30°: 


fiH 


TABLE A -2 


Curve Kit Conparinon 


F / 9012 P 

■^data' 'calculated/. 9012 


0 

1.0 

1.0 

4.4753 

0.9819 

0.9824 

8.7425 

0.9294 

0.9293 

12.5642 

0.8486 

0.8452 

15.6465 

0.7483 

0.7483 

20. 

— 

0.5651 

25. 

— 

0.2974 

30. 

— 

0. 


Fig. A-2 illustrates the zero sun cone angle case for positive 9 
and using the curve fit of Eq. (A.l) to extrapolate the data to e«30°. 
A curve is shown for the sun cone angle of 15° also. No fit of the 
data to a mathematical expression was attempted for that case. 



THRUST CONE ANGLE (ilegt 

Figure A-2 Force magnitude variation with cone angle 
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Table B-2 Partlals of and Yj^ with 
Respect to h and k 
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Non-zero Partlals of M wUli Respect to q 
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Table B-O Second ParUals of .and Y"! wUh Ileapect to h and k 
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APPBKDXX C . 

(( 

SINGIjE AVKHAGED OBLATENESS EQUATIONS 


ift these equations It, is the earth's radiua, |i is dte gravita- 
tional constantf n is the orbital angular speed, Jg the oblateness 
coefficient and a, h, k, p, q the equinoctial orbital elements. 


{[ 

Table 0-1 Variation of Parameters Equations | 
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Table C-2 Partial of Jg Equations with Respect to a 
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Table C-5 Partial of ISquations with Respect to p 
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Thitj appendix essentially repyoduces the summary of shadowing 
results from Uof. 12, From geometrical considerations an equation 
can bo derived which the entry and exit angles must satisfy. Such 
an equation is given in Reference 28 and the equation given in this 
section ia essentially the same; except that it is given in terms 
of equinoctial orbital elements. 

The spacecraft position is given by 

r = 

where X'j^ and were given in Eq. (3.16) and (3.17). hot the unit 
vector from the planet to the sun be given by 

Eg - + Zg£ 

This is in terms of the equinoctial coordinate frome and t!ms depends 
on the equinoctial orbital elements p and q. The calculation of 
the sun's direction in the equinoctial coordinate system is discussed 
in Section 3. If designates the planet's radius, the cosine of 
the angle between r and R is given by 


— s 


(|r|^ 


2 _ , 2 . 1/2 


Vs * " V 


Squaring and rearranging 


S 3 (l-x^jx^ H. (l-v2,v2 . 2X^V^X^Vj - = 0 





This is the shadow equation which must be satisfied by the entry 
and exit angles. and are functions of cosF, sinF, a, h, and k 
(see Eq. (3.16) and (3.17). By further manipulations one can derive 
a quartic equation in cosF. The coefficients of this quartic equation 
are given in Table D~l. Spurious roots can be eliminated by the 
criteria that S * 0 and that R * £ < 0» In addition, for the entry 
angle 9S/3F < 0 and for the exit angle 3S/8P >0. 


Table D-1 The Shadow Quartic Equation 
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Tha dorivaliivg of F v/ith tospect to x is needed to evaluate the 
coatate equation, lb can be obtained implicitly from tha ehadov; 
equation. 


Dm dp 


(D,5) 


These partials are listed in Table D-2 . Note that in calculatiny 
DS/Dp and 9s/Dq v/e have taken into account the fact that the sun's 
direction is given in equinoctial coordinates. 


Table b“2 Partlals nl' the Shadow Function 
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APPENDIX 13 



A a'l’UATEGY WHICH MAXIMIZES ^^EAM imTE-OP-CHANGE OF EKERGV- , 


,, Ai3 part oJ; a solar sail planetocentric performance analysis, a 
strategy was developed, cpded/ and run to produce trajectories which 
maximize the mean orbital rate-of-change of energy locally, The 
computer program {SUNSPOT} , which was being developed to produce mini- 
mum time escape trajectories, w’as modified to produce a number of 

I'.. 

trajectories using this strategy. 

The computer code, which uses the method of averaging, is valid 

for the spiral phase only. Tije averaging method is valid only if 

V ’’ 

the thrust/weight ratio is small. Therefore the trajectories cannot 
be extended to escape energy (C^ - 0) • 

Although the computer program is- designed to -solve a- tv.’o point 
boundary “value problem, single trajectories can be produced. The 
state and costate equations- are integrated, starting at a specified 
initial orbit, for a specified length of time. The state consists 
of -five averaged orbital elriments; the costate consissts of the ad- 
joints to these elements. In order to produce traj e,t tor ies using 
the maximum rate~of-change of energy strategy, the costate equations 
were not used. Instead the adjoint to the semimajor axis was set 
to 6.’ nonzero value, and the other ad joints were constrained to be 
zero, jlinde the eyiergy is proportional to the negative of the 
inverse of the semimajor axis, this strategy has the desired result. 
This strategy is equivaleht to forcing the primer vector to be 
tangent to the orbit at all times. (For the solar sail, this does 
not mean that the thrust is tangential) . 


The strategy requires that the energy change be maximized where 
the energy, E, is given by 


"U ra 
a 2 


(E.l) 


Thus 

jImNK CUMP 


35 


2 
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E e 


(E.2) 



tt 


where in ishha mass (consbant) , w is the gravitational constanti and 
a is the semimajor axis. Thus E is maximized by maximizing a. Since 
the Variation of par-ameters equation for semimajor axial is given by 

-i • 

^ ^ 2 S-. u (E.3) 

P 


Then a is maximized when u has the largest projection onto the 
•.‘elocity vector, V. For the sail, the thrust direction, u, will not 
necessarily be colinear with V, The above procedure is simply ac- 
complished using the optimization program by constraining the adjoint 
to semimajor axis to be a non-zero Constant and the other adjoints to 
be zero. The variation of parameters equations for all the orbital 
elements, z, can be written (M is a matrix): 


Then the Hamiltonian is 



(E,4) 


, T • 
H == >:_ Z 


(E.S) 


» i T 

Then if the ad;]oxnt vecbor is X 


[k, 0^}, H = [k, Z = k a 


Thus maximizing H, maximizes a. Since the optimizing program uses 
the method of averaging, it is the mean energy rate that is ac- 
tually maxiTTiized. 

The thrust vector for the solar sail is calculated as a function 

T 

of the primer vector (or equivalently, M M . For a flat perfectly- 

reflecting sail, the force acts normal to the sail and is propor- 

2 . 

tional to cos 0, where 6, the cone angle, is tne angle between the 
normal and the sun - vehicle vector. The relationship between 
thrust direction and primer vector has been previously derived, 
for example, in Ref. 6. A more accurate force model for the square 
sail is given Ref. 22. In this case the force is proportional to 
t C 2 cos 2 G -i- cos 43 , where the cone angle, S, is the angle 
between the force direction and the sun-vehicle vector. In this 
case the calculation of the thrust vector is m.'re involved, but 
has also been implemented in the SUNSPOT code. 


SB 


' A number ofc trajectories were generated using the niaxlnium 
ratu-of-change oi; energy strategy. Different initial orbits woro 
assumed; some cases used the idealiae( '•'il, others tiio more accurate 
model; soiiio casos included shadowing. case did not includo solar 

motion. Oblateness. was not included in any of the cases. A memo 
describing these cases was sent to the technical monitor . 
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